home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / gamma.i < prev    next >
Text File  |  1995-07-26  |  2KB  |  83 lines
  1. /*
  2.    GAMMA.I
  3.    Gamma function, beta function, and binomial coefficients.
  4.  
  5.    $Id$
  6.  */
  7. /*    Copyright (c) 1994.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. /* Algorithms from Numerical Recipes, Press, et. al., section 6.1,
  11.    Cambridge University Press (1988).
  12.  */
  13.  
  14. func ln_gamma(z)
  15. /* DOCUMENT ln_gamma(z)
  16.      returns natural log of the gamma function.
  17.      Error is less than 2.e-10 for real part of z>=1.
  18.      Use lngamma if you know that all z>=1, or if you don't care much about
  19.      accuracy for z<1.
  20.    SEE ALSO: lngamma, bico, beta
  21.  */
  22. {
  23.   if (structof(z)==complex) big= (z.re>=1.0);
  24.   else big= (z>=1.0);
  25.   list= where(big);
  26.   if (numberof(list)) lg1= lngamma(z(list));
  27.   list= where(!big);
  28.   if (numberof(list)) {
  29.     z= z(list);
  30.     lg2= log(pi/sin(pi*z)) - lngamma(1.0-z)
  31.   }
  32.   return merge(lg1,lg2,big);
  33. }
  34.  
  35. func lngamma(x)
  36. /* DOCUMENT lngamma(x)
  37.      returns natural log of the gamma function.
  38.      Error is less than 2.e-10 for real part of x>=1.
  39.      Use ln_gamma if some x<1.
  40.    SEE ALSO: ln_gamma, bico, beta
  41.  */
  42. {
  43.   tmp= x+4.5;
  44.   tmp-= (x-0.5)*log(tmp);
  45.   ser= 1.0 + 76.18009173/x - 86.50532033/(x+1.) + 24.01409822/(x+2.) -
  46.     1.231739516/(x+3.) + 0.120858003e-2/(x+4.) - 0.536382e-5/(x+5.);
  47.   return -tmp+log(2.50662827465*ser);
  48. }
  49.  
  50. func bico(n,k)
  51. /* DOCUMENT bico(n,k)
  52.      returns the binomial coefficient n!/(k!(n-k)!) as a double.
  53.    SEE ALSO: ln_gamma, beta
  54.  */
  55. {
  56.   return floor(0.5+exp(lngamma(n+1.)-lngamma(k+1.)-lngamma(n-k+1.)));
  57. }
  58.  
  59. func beta(z,w)
  60. /* DOCUMENT beta(z,w)
  61.      returns the beta function gamma(z)gamma(w)/gamma(z+w)
  62.    SEE ALSO: ln_gamma, bico
  63.  */
  64. {
  65.   return exp(lngamma(z)+lngamma(w)-lngamma(z+w));
  66. }
  67.  
  68. func erfc(x)
  69. /* DOCUMENT erfc(x)
  70.      returns the complementary error function 1-erf with fractional
  71.      error less than 1.2e-7 everywhere.
  72.  */
  73. {
  74.   z= abs(x);
  75.   t= 1.0/(1.0+0.5*z);
  76.   ans= t*exp(-z*z +
  77.          poly(t, -1.26551223, 1.00002368, 0.37409196, 0.09678418,
  78.           -0.18628806, 0.27886807, -1.13520398, 1.48851587,
  79.           -0.82215223, 0.17087277));
  80.   z= sign(x);
  81.   return ans*z + (1.-z);
  82. }
  83.